library(MASS)
library(gplots)
library(ggplot2)
library(Hmisc)
library(sciplot)
library(RItools)
library(aod)
library(memisc)
library(lme4)
library(arm)
library(gmodels)

remove(list=ls())

###
# Treatment Effects, Beta Figures
###

setwd("/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/Kstan_Data")
d <- read.csv(file="kstan_r_2016_05_21.csv",head=TRUE,sep=",")

d <- d[ which(d$Q15 != "Difficult to answer/Refuse"), ]



###	Beta Figure: Investigation DV, Full Sample (Excluding DKRTA)
###		Figure 2 in manuscript
nbyti<- as.matrix(aggregate(approve_inv ~ tmt , data=d, FUN=function(x) sum(!is.na(x)) ))
	n1<-nbyti[1,2]
	n2<-nbyti[2,2]
d2 <- d[ which(d$approve_inv =="1"), ]
abyti<- as.matrix(aggregate(approve_inv ~ tmt, data=d2, FUN=function(x) sum(!is.na(x)) ))
	a1<-abyti[1,2]
	a2<-abyti[2,2]
n<-5000
pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
tis <- as.matrix(rbind(t1,t2))
pis <- as.matrix(rbind(pi_1,pi_2))
betas <- as.data.frame(cbind(tis,pis))
betas$treatmentname = factor(betas$V1, labels = c("Control","Treatment"))
lineplot.CI(x.factor = betas$treatmentname, response = betas$V2, type="p",
	data = betas, main = "", xlab = "Treatment Group", ylab = "% Approving Investigation",
	ylim=c(0.60,0.9),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_nodk.eps", horizontal=FALSE)



###	Beta Figure: Investigation DV, Non-OshObJal Sample
###		Figure 3 in manuscript
d_nonosh <- d[ which(d$oshobjal=="0"), ]
nbyti_n<- as.matrix(aggregate(approve_inv ~ tmt , data=d_nonosh, FUN=function(x) sum(!is.na(x)) ))
	n1_n<-nbyti_n[1,2]
	n2_n<-nbyti_n[2,2]
d2_n <- d_nonosh[ which(d_nonosh$approve_inv =="1"), ]
abyti_n<- as.matrix(aggregate(approve_inv ~ tmt, data=d2_n, FUN=function(x) sum(!is.na(x)) ))
	a1_n<-abyti_n[1,2]
	a2_n<-abyti_n[2,2]
n_n<-5000
pi_1_n <- as.matrix(rbeta(n_n, a1_n+0.5, n1_n-a1_n+0.5))
pi_2_n <- as.matrix(rbeta(n_n, a2_n+0.5, n2_n-a2_n+0.5))
t1_n <- as.matrix(rep(1,n_n, nrow=n_n, ncol=1))
t2_n <- as.matrix(rep(2,n_n, nrow=n_n, ncol=1))
tis_n <- as.matrix(rbind(t1_n,t2_n))
pis_n <- as.matrix(rbind(pi_1_n,pi_2_n))
betas_n <- as.data.frame(cbind(tis_n,pis_n))
betas_n$treatmentname = factor(betas_n$V1, labels = c("Control","Treatment"))
lineplot.CI(x.factor = betas_n$treatmentname, response = betas_n$V2, type="p",
	data = betas_n, main = "", xlab = "Treatment Group", ylab = "% Approving Investigation",
	ylim=c(0.60,0.90),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_nonoshobjal_nodk.eps", horizontal=FALSE)


###	Beta Figure: Investigation DV, OshObJal Sample
###		Figure 3 in manuscript
d_osh <- d[ which(d$oshobjal=="1"), ]
nbyti<- as.matrix(aggregate(approve_inv ~ tmt , data=d_osh, FUN=function(x) sum(!is.na(x)) ))
	n1<-nbyti[1,2]
	n2<-nbyti[2,2]
d2 <- d_osh[ which(d_osh$approve_inv =="1"), ]
abyti<- as.matrix(aggregate(approve_inv ~ tmt, data=d2, FUN=function(x) sum(!is.na(x)) ))
	a1<-abyti[1,2]
	a2<-abyti[2,2]
n<-5000
pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
tis <- as.matrix(rbind(t1,t2))
pis <- as.matrix(rbind(pi_1,pi_2))
betas <- as.data.frame(cbind(tis,pis))
betas$treatmentname = factor(betas$V1, labels = c("Control","Treatment"))
lineplot.CI(x.factor = betas$treatmentname, response = betas$V2, type="p",
	data = betas, main = "", xlab = "Treatment Group", ylab = "% Approving Investigation",
	ylim=c(0.60,0.90),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_oshobjal_nodk.eps", horizontal=FALSE)


par(mfrow=c(1,2))

lineplot.CI(x.factor = betas_n$treatmentname, response = betas_n$V2, type="p",
	data = betas_n, main = "", xlab = "Non-Osh/Ob./Jal.", ylab = "% Approving Investigation",
	ylim=c(0.60,0.90),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
lineplot.CI(x.factor = betas$treatmentname, response = betas$V2, type="p",
	data = betas, main = "", xlab = "Osh/Ob./Jal.", ylab = "% Approving Investigation",
	ylim=c(0.60,0.90),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_combo_nodk.eps", horizontal=FALSE)



###
# Treatment effects, Uzbek/NonUzebek
#		Figures 5 and 6 in manuscript
###

remove(list=ls())
d <- read.csv(file="kstan_r_2016_05_21.csv",head=TRUE,sep=",")
	d <- d[ which(d$Q15 != "Difficult to answer/Refuse"), ]
	d_uzbek <- d[ which(d$uzbek_nat=="1"), ]
	d_nonuzbek <- d[ which(d$uzbek_nat=="0"), ]

###	Beta Figure: Investigation DV, Full Sample, Uzbek Only (Excluding DKRTA)
nbyti<- as.matrix(aggregate(approve_inv ~ tmt , data=d_uzbek, FUN=function(x) sum(!is.na(x)) ))
	n1<-nbyti[1,2]
	n2<-nbyti[2,2]
d2 <- d_uzbek[ which(d_uzbek$approve_inv =="1"), ]
abyti<- as.matrix(aggregate(approve_inv ~ tmt, data=d2, FUN=function(x) sum(!is.na(x)) ))
	a1<-abyti[1,2]
	a2<-abyti[2,2]
n<-5000
pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
tis <- as.matrix(rbind(t1,t2))
pis <- as.matrix(rbind(pi_1,pi_2))
betas_uzbek <- as.data.frame(cbind(tis,pis))
betas_uzbek$treatmentname = factor(betas_uzbek$V1, labels = c("Control","Treatment"))

lineplot.CI(x.factor = betas_uzbek$treatmentname, response = betas_uzbek$V2, type="p",
	data = betas_uzbek, main = "Uzbek, All Regions", xlab = "Treatment Group", ylab = "% Approving Investigation",
	ylim=c(0.50,1.0),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_nodk_uzbek_allregion.eps", horizontal=FALSE)

###	Beta Figure: Investigation DV, Full Sample, Non-Uzbek Only (Excluding DKRTA)
nbyti<- as.matrix(aggregate(approve_inv ~ tmt , data=d_nonuzbek , FUN=function(x) sum(!is.na(x)) ))
	n1<-nbyti[1,2]
	n2<-nbyti[2,2]
d2 <- d_nonuzbek [ which(d_nonuzbek $approve_inv =="1"), ]
abyti<- as.matrix(aggregate(approve_inv ~ tmt, data=d2, FUN=function(x) sum(!is.na(x)) ))
	a1<-abyti[1,2]
	a2<-abyti[2,2]
n<-5000
pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
tis <- as.matrix(rbind(t1,t2))
pis <- as.matrix(rbind(pi_1,pi_2))
betas_nonuzbek <- as.data.frame(cbind(tis,pis))
betas_nonuzbek $treatmentname = factor(betas_nonuzbek $V1, labels = c("Control","Treatment"))

lineplot.CI(x.factor = betas_nonuzbek $treatmentname, response = betas_nonuzbek $V2, type="p",
	data = betas_nonuzbek , main = "Non-Uzbek, All Regions", xlab = "Treatment Group", ylab = "% Approving Investigation",
	ylim=c(0.50,1.0),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_nodk_nonuzbek_allregion.eps", horizontal=FALSE)


par(mfrow=c(1,2))

lineplot.CI(x.factor = betas_uzbek$treatmentname, response = betas_uzbek$V2, type="p",
	data = betas_uzbek, main = "Uzbek, All Regions", xlab = "", ylab = "% Approving Investigation",
	ylim=c(0.40,1.0),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
lineplot.CI(x.factor = betas_nonuzbek $treatmentname, response = betas_nonuzbek $V2, type="p",
	data = betas_nonuzbek , main = "Non-Uzbek, All Regions", xlab = "", ylab = "% Approving Investigation",
	ylim=c(0.40,1.0),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_uzcombo_fullregion_nodk.eps", horizontal=FALSE)


###
# Treatment effects, Uzbek/NonUzebek, OOJ Only
###

	d_osh <- d[ which(d$oshobjal=="1"), ]
	d_uzbek <- d_osh[ which(d_osh$uzbek_nat=="1"), ]
	d_nonuzbek <- d_osh[ which(d_osh$uzbek_nat=="0"), ]

###	Beta Figure: Investigation DV, Full Sample, Uzbek Only (Excluding DKRTA)
nbyti<- as.matrix(aggregate(approve_inv ~ tmt , data=d_uzbek, FUN=function(x) sum(!is.na(x)) ))
	n1<-nbyti[1,2]
	n2<-nbyti[2,2]
d2 <- d_uzbek[ which(d_uzbek$approve_inv =="1"), ]
abyti<- as.matrix(aggregate(approve_inv ~ tmt, data=d2, FUN=function(x) sum(!is.na(x)) ))
	a1<-abyti[1,2]
	a2<-abyti[2,2]
n<-5000
pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
tis <- as.matrix(rbind(t1,t2))
pis <- as.matrix(rbind(pi_1,pi_2))
betas_uzbek <- as.data.frame(cbind(tis,pis))
betas_uzbek$treatmentname = factor(betas_uzbek$V1, labels = c("Control","Treatment"))

lineplot.CI(x.factor = betas_uzbek$treatmentname, response = betas_uzbek$V2, type="p",
	data = betas_uzbek, main = "Uzbek, Osh/Osh Oblast/Jalal-Abad", xlab = "Treatment Group", ylab = "% Approving Investigation",
	ylim=c(0.50,1.0),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_nodk_uzbek_ooj.eps", horizontal=FALSE)

###	Beta Figure: Investigation DV, Full Sample, Non-Uzbek Only (Excluding DKRTA)
nbyti<- as.matrix(aggregate(approve_inv ~ tmt , data=d_nonuzbek , FUN=function(x) sum(!is.na(x)) ))
	n1<-nbyti[1,2]
	n2<-nbyti[2,2]
d2 <- d_nonuzbek [ which(d_nonuzbek $approve_inv =="1"), ]
abyti<- as.matrix(aggregate(approve_inv ~ tmt, data=d2, FUN=function(x) sum(!is.na(x)) ))
	a1<-abyti[1,2]
	a2<-abyti[2,2]
n<-5000
pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
tis <- as.matrix(rbind(t1,t2))
pis <- as.matrix(rbind(pi_1,pi_2))
betas_nonuzbek <- as.data.frame(cbind(tis,pis))
betas_nonuzbek $treatmentname = factor(betas_nonuzbek $V1, labels = c("Control","Treatment"))

lineplot.CI(x.factor = betas_nonuzbek $treatmentname, response = betas_nonuzbek $V2, type="p",
	data = betas_nonuzbek , main = "Non-Uzbek, Osh/Osh Oblast/Jalal-Abad", xlab = "Treatment Group", ylab = "% Approving Investigation",
	ylim=c(0.50,1.0),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_nodk_nonuzbek_ooj.eps", horizontal=FALSE)


par(mfrow=c(1,2))

lineplot.CI(x.factor = betas_uzbek$treatmentname, response = betas_uzbek$V2, type="p",
	data = betas_uzbek, main = "Uzbek", xlab = "", ylab = "% Approving Investigation",
	ylim=c(0.40,1.0),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
lineplot.CI(x.factor = betas_nonuzbek $treatmentname, response = betas_nonuzbek $V2, type="p",
	data = betas_nonuzbek , main = "Non-Uzbek", xlab = "", ylab = "% Approving Investigation",
	ylim=c(0.40,1.0),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_uzcombo_ooj_nodk.eps", horizontal=FALSE)


###
# Treatment effects, Uzbek/NonUzebek, non-OOJ Only, R3 request
###

	d_notosh <- d[ which(d$oshobjal=="0"), ]
	d_uzbek <- d_notosh[ which(d_notosh$uzbek_nat=="1"), ]
	d_nonuzbek <- d_notosh[ which(d_notosh$uzbek_nat=="0"), ]

###	Beta Figure: Investigation DV, Full Sample, Uzbek Only (Excluding DKRTA)
nbyti<- as.matrix(aggregate(approve_inv ~ tmt , data=d_uzbek, FUN=function(x) sum(!is.na(x)) ))
	n1<-nbyti[1,2]
	n2<-nbyti[2,2]
d2 <- d_uzbek[ which(d_uzbek$approve_inv =="1"), ]
abyti<- as.matrix(aggregate(approve_inv ~ tmt, data=d2, FUN=function(x) sum(!is.na(x)) ))
	a1<-abyti[1,2]
	a2<-abyti[2,2]
n<-5000
pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
tis <- as.matrix(rbind(t1,t2))
pis <- as.matrix(rbind(pi_1,pi_2))
betas_uzbek <- as.data.frame(cbind(tis,pis))
betas_uzbek$treatmentname = factor(betas_uzbek$V1, labels = c("Control","Treatment"))

lineplot.CI(x.factor = betas_uzbek$treatmentname, response = betas_uzbek$V2, type="p",
	data = betas_uzbek, main = "Uzbek, non-Osh", xlab = "Treatment Group", ylab = "% Approving Investigation",
	ylim=c(0.50,1.0),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_nodk_uzbek_nonooj.eps", horizontal=FALSE)

###	Beta Figure: Investigation DV, Full Sample, Non-Uzbek Only (Excluding DKRTA)
nbyti<- as.matrix(aggregate(approve_inv ~ tmt , data=d_nonuzbek , FUN=function(x) sum(!is.na(x)) ))
	n1<-nbyti[1,2]
	n2<-nbyti[2,2]
d2 <- d_nonuzbek [ which(d_nonuzbek $approve_inv =="1"), ]
abyti<- as.matrix(aggregate(approve_inv ~ tmt, data=d2, FUN=function(x) sum(!is.na(x)) ))
	a1<-abyti[1,2]
	a2<-abyti[2,2]
n<-5000
pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
tis <- as.matrix(rbind(t1,t2))
pis <- as.matrix(rbind(pi_1,pi_2))
betas_nonuzbek <- as.data.frame(cbind(tis,pis))
betas_nonuzbek $treatmentname = factor(betas_nonuzbek $V1, labels = c("Control","Treatment"))

lineplot.CI(x.factor = betas_nonuzbek $treatmentname, response = betas_nonuzbek $V2, type="p",
	data = betas_nonuzbek , main = "Non-Uzbek, non-Osh", xlab = "Treatment Group", ylab = "% Approving Investigation",
	ylim=c(0.50,1.0),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_nodk_nonuzbek_nonooj.eps", horizontal=FALSE)


par(mfrow=c(1,2))

lineplot.CI(x.factor = betas_uzbek$treatmentname, response = betas_uzbek$V2, type="p",
	data = betas_uzbek, main = "Uzbek", xlab = "", ylab = "% Approving Investigation",
	ylim=c(0.40,1.0),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
lineplot.CI(x.factor = betas_nonuzbek $treatmentname, response = betas_nonuzbek $V2, type="p",
	data = betas_nonuzbek , main = "Non-Uzbek", xlab = "", ylab = "% Approving Investigation",
	ylim=c(0.40,1.0),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_uzcombo_nonooj_nodk.eps", horizontal=FALSE)





###
# Treatment effects, Gov App/Non Gov App
#		Figure 12 in Appendix
###

# data$approve_gov1 : 1,2 (do not approve) 3,4 (approve)

remove(list=ls())
d <- read.csv(file="kstan_r_2016_05_21.csv",head=TRUE,sep=",")
	d <- d[ which(d$Q15 != "Difficult to answer/Refuse"), ]
	d_govnoapp <- d[ which(d$approve_gov1 =="1" | d$approve_gov1 =="2"), ]
	d_govapp <- d[ which(d$approve_gov1 =="3" | d$approve_gov1 =="4"), ]

###	Beta Figure: Investigation DV, Full Sample, Gov App Only (Excluding DKRTA)
nbyti<- as.matrix(aggregate(approve_inv ~ tmt , data= d_govapp, FUN=function(x) sum(!is.na(x)) ))
	n1<-nbyti[1,2]
	n2<-nbyti[2,2]
d2 <- d_govapp[ which(d_govapp $approve_inv =="1"), ]
abyti<- as.matrix(aggregate(approve_inv ~ tmt, data=d2, FUN=function(x) sum(!is.na(x)) ))
	a1<-abyti[1,2]
	a2<-abyti[2,2]
n<-5000
pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
tis <- as.matrix(rbind(t1,t2))
pis <- as.matrix(rbind(pi_1,pi_2))
betas_govapp <- as.data.frame(cbind(tis,pis))
betas_govapp$treatmentname = factor(betas_govapp$V1, labels = c("Control","Treatment"))

lineplot.CI(x.factor = betas_govapp$treatmentname, response = betas_govapp$V2, type="p",
	data = betas_govapp, main = "Approve of Gov.", xlab = "Treatment Group", ylab = "% Approving Investigation",
	ylim=c(0.50,1.0),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_nodk_govapp_allregion.eps", horizontal=FALSE)

###	Beta Figure: Investigation DV, Full Sample, No Gov App Only (Excluding DKRTA)
nbyti<- as.matrix(aggregate(approve_inv ~ tmt , data= d_govnoapp , FUN=function(x) sum(!is.na(x)) ))
	n1<-nbyti[1,2]
	n2<-nbyti[2,2]
d2 <- d_govnoapp [ which(d_govnoapp$approve_inv =="1"), ]
abyti<- as.matrix(aggregate(approve_inv ~ tmt, data=d2, FUN=function(x) sum(!is.na(x)) ))
	a1<-abyti[1,2]
	a2<-abyti[2,2]
n<-5000
pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
tis <- as.matrix(rbind(t1,t2))
pis <- as.matrix(rbind(pi_1,pi_2))
betas_govnoapp <- as.data.frame(cbind(tis,pis))
betas_govnoapp$treatmentname = factor(betas_govnoapp$V1, labels = c("Control","Treatment"))

lineplot.CI(x.factor = betas_govnoapp$treatmentname, response = betas_govnoapp$V2, type="p",
	data = betas_govnoapp , main = "Not Approve of Gov.", xlab = "Treatment Group", ylab = "% Approving Investigation",
	ylim=c(0.50,1.0),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_nodk_nogovapp_allregion.eps", horizontal=FALSE)


par(mfrow=c(1,2))

lineplot.CI(x.factor = betas_govapp$treatmentname, response = betas_govapp$V2, type="p",
	data = betas_govapp, main = "Approve of Gov.", xlab = "", ylab = "% Approving Investigation",
	ylim=c(0.40,1.0),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
lineplot.CI(x.factor = betas_govnoapp$treatmentname, response = betas_govnoapp$V2, type="p",
	data = betas_govnoapp , main = "Not Approve of Gov.", xlab = "", ylab = "% Approving Investigation",
	ylim=c(0.40,1.0),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_govappcombo_fullregion_nodk.eps", horizontal=FALSE)

tmgovapp <- betas_govapp[ which(betas_govapp$treatmentname == "Treatment"), ]
summary(tmgovapp$V2)
ctgovapp <- betas_govapp[ which(betas_govapp$treatmentname == "Control"), ]
summary(ctgovapp$V2)

tmgovnoapp <- betas_govnoapp[ which(betas_govnoapp$treatmentname == "Treatment"), ]
summary(tmgovnoapp$V2)
ctgovnoapp <- betas_govnoapp[ which(betas_govnoapp$treatmentname == "Control"), ]
summary(ctgovnoapp$V2)




###
# Treatment effects, Heard of ICC/Not Heard
#		Figure 7 in Manuscript
###

remove(list=ls())
d <- read.csv(file="kstan_r_2016_05_21.csv",head=TRUE,sep=",")
	d <- d[ which(d$Q15 != "Difficult to answer/Refuse"), ]
	d_heard <- d[ which(d$heard_icc=="1"), ]
	d_noheard <- d[ which(d$heard_icc=="0"), ]

###	Beta Figure: Investigation DV, Full Sample, Heard Only (Excluding DKRTA)
nbyti<- as.matrix(aggregate(approve_inv ~ tmt , data= d_heard, FUN=function(x) sum(!is.na(x)) ))
	n1<-nbyti[1,2]
	n2<-nbyti[2,2]
d2 <- d_heard[ which(d_heard$approve_inv =="1"), ]
abyti<- as.matrix(aggregate(approve_inv ~ tmt, data=d2, FUN=function(x) sum(!is.na(x)) ))
	a1<-abyti[1,2]
	a2<-abyti[2,2]
n<-5000
pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
tis <- as.matrix(rbind(t1,t2))
pis <- as.matrix(rbind(pi_1,pi_2))
betas_heard <- as.data.frame(cbind(tis,pis))
betas_heard$treatmentname = factor(betas_heard$V1, labels = c("Control","Treatment"))

lineplot.CI(x.factor = betas_heard$treatmentname, response = betas_heard$V2, type="p",
	data = betas_heard, main = "Heard of ICC", xlab = "Treatment Group", ylab = "% Approving Investigation",
	ylim=c(0.50,1.0),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_nodk_heard_allregion.eps", horizontal=FALSE)




###	Beta Figure: Investigation DV, Full Sample, No Heard Only (Excluding DKRTA)
nbyti<- as.matrix(aggregate(approve_inv ~ tmt , data= d_noheard , FUN=function(x) sum(!is.na(x)) ))
	n1<-nbyti[1,2]
	n2<-nbyti[2,2]
d2 <- d_noheard [ which(d_noheard$approve_inv =="1"), ]
abyti<- as.matrix(aggregate(approve_inv ~ tmt, data=d2, FUN=function(x) sum(!is.na(x)) ))
	a1<-abyti[1,2]
	a2<-abyti[2,2]
n<-5000
pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
tis <- as.matrix(rbind(t1,t2))
pis <- as.matrix(rbind(pi_1,pi_2))
betas_noheard <- as.data.frame(cbind(tis,pis))
betas_noheard$treatmentname = factor(betas_noheard$V1, labels = c("Control","Treatment"))

lineplot.CI(x.factor = betas_noheard$treatmentname, response = betas_noheard$V2, type="p",
	data = betas_noheard , main = "Not Heard of ICC", xlab = "Treatment Group", ylab = "% Approving Investigation",
	ylim=c(0.50,1.0),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_nodk_noheard_allregion.eps", horizontal=FALSE)


par(mfrow=c(1,2))

lineplot.CI(x.factor = betas_heard$treatmentname, response = betas_heard$V2, type="p",
	data = betas_heard, main = "Heard of ICC", xlab = "", ylab = "% Approving Investigation",
	ylim=c(0.40,1.0),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
lineplot.CI(x.factor = betas_noheard$treatmentname, response = betas_noheard$V2, type="p",
	data = betas_noheard , main = "Not Heard of ICC", xlab = "", ylab = "% Approving Investigation",
	ylim=c(0.40,1.0),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_heardcombo_fullregion_nodk.eps", horizontal=FALSE)

tmheard <- betas_heard[ which(betas_heard$treatmentname == "Treatment"), ]
summary(tmheard $V2)
ctheard <- betas_heard[ which(betas_heard$treatmentname == "Control"), ]
summary(ctgovapp$V2)

tmnoheard <- betas_noheard[ which(betas_noheard$treatmentname == "Treatment"), ]
summary(tmnoheard $V2)
ctnoheard <- betas_govnoapp[ which(betas_noheard$treatmentname == "Control"), ]
summary(ctnoheard$V2)







###
# Treatment Effects, Beta Figures, ICC DV, Appendix Fig 7
###

remove(list=ls())
setwd("/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/Kstan_Data")
d <- read.csv(file="kstan_r_2016_05_21.csv",head=TRUE,sep=",")
	# Also need to change the data call for the balance tests at bottom
d <- d[ which(d$Q16 != "Difficult to answer/Refuse"), ]


###	Beta Figure: ICC DV, Full Sample
nbyti<- as.matrix(aggregate(approve_icc ~ tmt , data=d, FUN=function(x) sum(!is.na(x)) ))
	n1<-nbyti[1,2]
	n2<-nbyti[2,2]
d2 <- d[ which(d$approve_icc =="1"), ]
abyti<- as.matrix(aggregate(approve_icc ~ tmt, data=d2, FUN=function(x) sum(!is.na(x)) ))
	a1<-abyti[1,2]
	a2<-abyti[2,2]
n<-5000
pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
tis <- as.matrix(rbind(t1,t2))
pis <- as.matrix(rbind(pi_1,pi_2))
betas <- as.data.frame(cbind(tis,pis))
betas$treatmentname = factor(betas$V1, labels = c("Control","Treatment"))
lineplot.CI(x.factor = betas$treatmentname, response = betas$V2, type="p",
	data = betas, main = "", xlab = "Treatment Group", ylab = "% Approving ICC",
	ylim=c(0.6,0.90),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveicc_beta_nodk.eps", horizontal=FALSE)


###	Beta Figure: ICC DV, Non-Osh Sample
d_nonosh <- d[ which(d$osh=="0"), ]
nbyti<- as.matrix(aggregate(approve_icc ~ tmt , data=d_nonosh, FUN=function(x) sum(!is.na(x)) ))
	n1<-nbyti[1,2]
	n2<-nbyti[2,2]
d2 <- d_nonosh[ which(d_nonosh$approve_icc =="1"), ]
abyti<- as.matrix(aggregate(approve_icc ~ tmt, data=d2, FUN=function(x) sum(!is.na(x)) ))
	a1<-abyti[1,2]
	a2<-abyti[2,2]
n<-5000
pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
tis <- as.matrix(rbind(t1,t2))
pis <- as.matrix(rbind(pi_1,pi_2))
betas <- as.data.frame(cbind(tis,pis))
betas$treatmentname = factor(betas$V1, labels = c("Control","NIMBY"))
lineplot.CI(x.factor = betas$treatmentname, response = betas$V2, type="p",
	data = betas, main = "", xlab = "Treatment Group", ylab = "% Approving ICC",
	ylim=c(0.40,0.55),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveicc_beta_nonosh.eps", horizontal=FALSE)


###	Beta Figure: ICC DV, Osh Sample
d_osh <- d[ which(d$osh=="1"), ]
nbyti<- as.matrix(aggregate(approve_icc ~ tmt , data=d_osh, FUN=function(x) sum(!is.na(x)) ))
	n1<-nbyti[1,2]
	n2<-nbyti[2,2]
d2 <- d_osh[ which(d_osh$approve_icc =="1"), ]
abyti<- as.matrix(aggregate(approve_icc ~ tmt, data=d2, FUN=function(x) sum(!is.na(x)) ))
	a1<-abyti[1,2]
	a2<-abyti[2,2]
n<-5000
pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
tis <- as.matrix(rbind(t1,t2))
pis <- as.matrix(rbind(pi_1,pi_2))
betas <- as.data.frame(cbind(tis,pis))
betas$treatmentname = factor(betas$V1, labels = c("Control","NIMBY"))
lineplot.CI(x.factor = betas$treatmentname, response = betas$V2, type="p",
	data = betas, main = "", xlab = "Treatment Group", ylab = "% ICC Investigation",
	ylim=c(0.55,0.85),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveicc_beta_osh.eps", horizontal=FALSE)




###
# WITH DKRTA, Appendix Figures 8-11
###

###	Beta Figures: Inv DV WITH DKRTA
remove(list=ls())
setwd("/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/Kstan_Data")
d <- read.csv(file="kstan_r_2016_05_21.csv",head=TRUE,sep=",")
	#approve_inv is fine for this, so long as you don't remove the observations with Q15 == DKRTA
	#d <- d[ which(d$Q15 != "Difficult to answer/Refuse"), ]

# Main treatment effect figure, WITH DKRTA
nbyti<- as.matrix(aggregate(approve_inv ~ tmt , data=d, FUN=function(x) sum(!is.na(x)) ))
	n1<-nbyti[1,2]
	n2<-nbyti[2,2]
d2 <- d[ which(d$approve_inv =="1"), ]
abyti<- as.matrix(aggregate(approve_inv ~ tmt, data=d2, FUN=function(x) sum(!is.na(x)) ))
	a1<-abyti[1,2]
	a2<-abyti[2,2]
n<-5000
pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
tis <- as.matrix(rbind(t1,t2))
pis <- as.matrix(rbind(pi_1,pi_2))
betas <- as.data.frame(cbind(tis,pis))
betas$treatmentname = factor(betas$V1, labels = c("Control","Treatment"))
lineplot.CI(x.factor = betas$treatmentname, response = betas$V2, type="p",
	data = betas, main = "", xlab = "Treatment Group", ylab = "% Approving Investigation",
	ylim=c(0.50,0.7),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_yesdk.eps", horizontal=FALSE)


###	Beta Figure: Investigation DV, Non-OshObJal Sample
d_nonosh <- d[ which(d$oshobjal=="0"), ]
nbyti_n<- as.matrix(aggregate(approve_inv ~ tmt , data=d_nonosh, FUN=function(x) sum(!is.na(x)) ))
	n1_n<-nbyti_n[1,2]
	n2_n<-nbyti_n[2,2]
d2_n <- d_nonosh[ which(d_nonosh$approve_inv =="1"), ]
abyti_n<- as.matrix(aggregate(approve_inv ~ tmt, data=d2_n, FUN=function(x) sum(!is.na(x)) ))
	a1_n<-abyti_n[1,2]
	a2_n<-abyti_n[2,2]
n_n<-5000
pi_1_n <- as.matrix(rbeta(n_n, a1_n+0.5, n1_n-a1_n+0.5))
pi_2_n <- as.matrix(rbeta(n_n, a2_n+0.5, n2_n-a2_n+0.5))
t1_n <- as.matrix(rep(1,n_n, nrow=n_n, ncol=1))
t2_n <- as.matrix(rep(2,n_n, nrow=n_n, ncol=1))
tis_n <- as.matrix(rbind(t1_n,t2_n))
pis_n <- as.matrix(rbind(pi_1_n,pi_2_n))
betas_n <- as.data.frame(cbind(tis_n,pis_n))
betas_n$treatmentname = factor(betas_n$V1, labels = c("Control","Treatment"))
lineplot.CI(x.factor = betas_n$treatmentname, response = betas_n$V2, type="p",
	data = betas_n, main = "", xlab = "Treatment Group", ylab = "% Approving Investigation",
	ylim=c(0.50,0.80),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_nonoshobjal_ydk.eps", horizontal=FALSE)


###	Beta Figure: Investigation DV, OshObJal Sample
d_osh <- d[ which(d$oshobjal=="1"), ]
nbyti<- as.matrix(aggregate(approve_inv ~ tmt , data=d_osh, FUN=function(x) sum(!is.na(x)) ))
	n1<-nbyti[1,2]
	n2<-nbyti[2,2]
d2 <- d_osh[ which(d_osh$approve_inv =="1"), ]
abyti<- as.matrix(aggregate(approve_inv ~ tmt, data=d2, FUN=function(x) sum(!is.na(x)) ))
	a1<-abyti[1,2]
	a2<-abyti[2,2]
n<-5000
pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
tis <- as.matrix(rbind(t1,t2))
pis <- as.matrix(rbind(pi_1,pi_2))
betas <- as.data.frame(cbind(tis,pis))
betas$treatmentname = factor(betas$V1, labels = c("Control","Treatment"))
lineplot.CI(x.factor = betas$treatmentname, response = betas$V2, type="p",
	data = betas, main = "", xlab = "Treatment Group", ylab = "% Approving Investigation",
	ylim=c(0.50,0.80),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_oshobjal_ydk.eps", horizontal=FALSE)

par(mfrow=c(1,2))
lineplot.CI(x.factor = betas_n$treatmentname, response = betas_n$V2, type="p",
	data = betas_n, main = "", xlab = "Non-Osh/Ob./Jal.", ylab = "% Approving Investigation",
	ylim=c(0.50,0.80),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
lineplot.CI(x.factor = betas$treatmentname, response = betas$V2, type="p",
	data = betas, main = "", xlab = "Osh/Ob./Jal.", ylab = "% Approving Investigation",
	ylim=c(0.50,0.80),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_combo_ydk.eps", horizontal=FALSE)



# Beta figures, treatment effects, Uzbek/NonUzebek, all regions, WITH DKRTA
remove(list=ls())
d <- read.csv(file="kstan_r_2016_05_21.csv",head=TRUE,sep=",")
	d_uzbek <- d[ which(d$uzbek_nat=="1"), ]
	d_nonuzbek <- d[ which(d$uzbek_nat=="0"), ]

###	Beta Figure: Investigation DV, Full Sample, Uzbek Only
nbyti<- as.matrix(aggregate(approve_inv ~ tmt , data=d_uzbek, FUN=function(x) sum(!is.na(x)) ))
	n1<-nbyti[1,2]
	n2<-nbyti[2,2]
d2 <- d_uzbek[ which(d_uzbek$approve_inv =="1"), ]
abyti<- as.matrix(aggregate(approve_inv ~ tmt, data=d2, FUN=function(x) sum(!is.na(x)) ))
	a1<-abyti[1,2]
	a2<-abyti[2,2]
n<-5000
pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
tis <- as.matrix(rbind(t1,t2))
pis <- as.matrix(rbind(pi_1,pi_2))
betas_uzbek <- as.data.frame(cbind(tis,pis))
betas_uzbek$treatmentname = factor(betas_uzbek$V1, labels = c("Control","Treatment"))

lineplot.CI(x.factor = betas_uzbek$treatmentname, response = betas_uzbek$V2, type="p",
	data = betas_uzbek, main = "Uzbek, All Regions", xlab = "Treatment Group", ylab = "% Approving Investigation",
	ylim=c(0.50,1.0),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_yesdk_uzbek_allregion.eps", horizontal=FALSE)

###	Beta Figure: Investigation DV, Full Sample, Non-Uzbek Only
nbyti<- as.matrix(aggregate(approve_inv ~ tmt , data=d_nonuzbek , FUN=function(x) sum(!is.na(x)) ))
	n1<-nbyti[1,2]
	n2<-nbyti[2,2]
d2 <- d_nonuzbek [ which(d_nonuzbek $approve_inv =="1"), ]
abyti<- as.matrix(aggregate(approve_inv ~ tmt, data=d2, FUN=function(x) sum(!is.na(x)) ))
	a1<-abyti[1,2]
	a2<-abyti[2,2]
n<-5000
pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
tis <- as.matrix(rbind(t1,t2))
pis <- as.matrix(rbind(pi_1,pi_2))
betas_nonuzbek <- as.data.frame(cbind(tis,pis))
betas_nonuzbek $treatmentname = factor(betas_nonuzbek $V1, labels = c("Control","Treatment"))

lineplot.CI(x.factor = betas_nonuzbek $treatmentname, response = betas_nonuzbek $V2, type="p",
	data = betas_nonuzbek , main = "Non-Uzbek, All Regions", xlab = "Treatment Group", ylab = "% Approving Investigation",
	ylim=c(0.50,1.0),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_yesdk_nonuzbek_allregion.eps", horizontal=FALSE)


par(mfrow=c(1,2))

lineplot.CI(x.factor = betas_uzbek$treatmentname, response = betas_uzbek$V2, type="p",
	data = betas_uzbek, main = "Uzbek, All Regions", xlab = "", ylab = "% Approving Investigation",
	ylim=c(0.30,0.8),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
lineplot.CI(x.factor = betas_nonuzbek $treatmentname, response = betas_nonuzbek $V2, type="p",
	data = betas_nonuzbek , main = "Non-Uzbek, All Regions", xlab = "", ylab = "% Approving Investigation",
	ylim=c(0.30,0.8),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_uzcombo_fullregion_yesdk.eps", horizontal=FALSE)


###
# Treatment effects, Uzbek/NonUzebek, OOJ Only
###

	d_osh <- d[ which(d$oshobjal=="1"), ]
	d_uzbek <- d_osh[ which(d_osh$uzbek_nat=="1"), ]
	d_nonuzbek <- d_osh[ which(d_osh$uzbek_nat=="0"), ]

###	Beta Figure: Investigation DV, Full Sample, Uzbek Only
nbyti<- as.matrix(aggregate(approve_inv ~ tmt , data=d_uzbek, FUN=function(x) sum(!is.na(x)) ))
	n1<-nbyti[1,2]
	n2<-nbyti[2,2]
d2 <- d_uzbek[ which(d_uzbek$approve_inv =="1"), ]
abyti<- as.matrix(aggregate(approve_inv ~ tmt, data=d2, FUN=function(x) sum(!is.na(x)) ))
	a1<-abyti[1,2]
	a2<-abyti[2,2]
n<-5000
pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
tis <- as.matrix(rbind(t1,t2))
pis <- as.matrix(rbind(pi_1,pi_2))
betas_uzbek <- as.data.frame(cbind(tis,pis))
betas_uzbek$treatmentname = factor(betas_uzbek$V1, labels = c("Control","Treatment"))

lineplot.CI(x.factor = betas_uzbek$treatmentname, response = betas_uzbek$V2, type="p",
	data = betas_uzbek, main = "Uzbek, Osh/Osh Oblast/Jalal-Abad", xlab = "Treatment Group", ylab = "% Approving Investigation",
	ylim=c(0.50,1.0),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_yesdk_uzbek_ooj.eps", horizontal=FALSE)

###	Beta Figure: Investigation DV, Full Sample, Non-Uzbek Only 
nbyti<- as.matrix(aggregate(approve_inv ~ tmt , data=d_nonuzbek , FUN=function(x) sum(!is.na(x)) ))
	n1<-nbyti[1,2]
	n2<-nbyti[2,2]
d2 <- d_nonuzbek [ which(d_nonuzbek $approve_inv =="1"), ]
abyti<- as.matrix(aggregate(approve_inv ~ tmt, data=d2, FUN=function(x) sum(!is.na(x)) ))
	a1<-abyti[1,2]
	a2<-abyti[2,2]
n<-5000
pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
tis <- as.matrix(rbind(t1,t2))
pis <- as.matrix(rbind(pi_1,pi_2))
betas_nonuzbek <- as.data.frame(cbind(tis,pis))
betas_nonuzbek $treatmentname = factor(betas_nonuzbek $V1, labels = c("Control","Treatment"))

lineplot.CI(x.factor = betas_nonuzbek $treatmentname, response = betas_nonuzbek $V2, type="p",
	data = betas_nonuzbek , main = "Non-Uzbek, Osh/Osh Oblast/Jalal-Abad", xlab = "Treatment Group", ylab = "% Approving Investigation",
	ylim=c(0.50,1.0),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_yesdk_nonuzbek_ooj.eps", horizontal=FALSE)


par(mfrow=c(1,2))

lineplot.CI(x.factor = betas_uzbek$treatmentname, response = betas_uzbek$V2, type="p",
	data = betas_uzbek, main = "Uzbek", xlab = "", ylab = "% Approving Investigation",
	ylim=c(0.30,0.8),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
lineplot.CI(x.factor = betas_nonuzbek $treatmentname, response = betas_nonuzbek $V2, type="p",
	data = betas_nonuzbek , main = "Non-Uzbek", xlab = "", ylab = "% Approving Investigation",
	ylim=c(0.30,0.8),
	ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_Kstan_Drafts/approveinv_beta_uzcombo_ooj_yesdk.eps", horizontal=FALSE)





###
# Bowers and Hansen Balance Tests, tables in Appendix A1.9, Tables 15 and 16
###
remove(list=ls())
setwd("/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/Kstan_Data")
d <- read.csv(file="kstan_r_2016_05_21.csv",head=TRUE,sep=",")

xBalance(tmt~ uzbek_nat + age_under50 + male + postsec_educ + employed + income_av, data=d,
	report="all")
	### Male is significant, and the overall statistic is significant

xBalance(tmt~ uzbek_nat + age_under50 + postsec_educ + employed + income_av, data=d,
	report="all")
	### Without male, the overall statistic is insignificant

xBalance(tmt~ region_1 + region_2 + region_3 + region_4 + region_5 + region_6 + region_7 + region_8 + region_9, data=d,
	report="all")
	### No region is significant, and the overall statistic is insignificant





d2 <- d[ which(d$Q15 != "Difficult to answer/Refuse"), ]
xBalance(tmt~ uzbek_nat + age_under50 + male + postsec_educ + employed + income_av, data=d2,
	report="all")
	### Male is significant, and the overall statistic is significant

xBalance(tmt~ uzbek_nat + age_under50 + postsec_educ + employed + income_av, data=d2,
	report="all")
	### Without male, the overall statistic is insignificant

xBalance(tmt~ region_1 + region_2 + region_3 + region_4 + region_5 + region_6 + region_7 + region_8 + region_9, data=d2,
	report="all")
	### No region is significant, and the overall statistic is insignificant








###
# MLM, Figure 4 in the manuscript
###
# copied and pasted from kstan_r_analysis_MLM_2016_05_25 on 3/18/2017


rm(list=ls())
setwd("/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/Kstan_Data")
data <- read.csv(file="kstan_r_2016_05_21.csv",head=TRUE,sep=",")

data <- data[ which(data$Q15 != "Difficult to answer/Refuse"), ]


n <- length(data$approve_inv)
y <- data$approve_inv
y_icc <- data$approve_icc
x <- data$tmt

region <- data$I1

### Multilevel Model: Treatment effects vary by region, with control variables

approve_gov1 <- data$approve_gov1
kyrgyzlang <- data$kyrgyzlang 
age_under50 <-data$age_under50
male  <-data$male
postsec_educ <-data$postsec_educ
employed  <-data$employed 
income_av <-data$income_av

# old specification, before we had approve as a moderator
#m2region <- glmer(y ~ x + (1+x|region) + approve_gov1 + kyrgyzlang + age_under50 + male + postsec_educ + employed + income_av, family = binomial)

m2region <- glmer(y ~ x + (1+x|region) + age_under50 + male + postsec_educ + employed + income_av, family = binomial)


m2region
ranef(m2region)
se.ranef(m2region)

model2<-mtable(m2region)
toLatex(model2)

### Figure XX: Region specific treatment effects

b.hat.m2 <- coef(m2region)$region[,2]
b.se.m2 <- se.coef(m2region)$region[,2]
J <- 9
regionvec <- as.numeric(1:J)

# Make vector of country names
region.name <- as.vector(data$I1)
uniq.name <- sort(unique(region.name))

lower <- as.matrix(b.hat.m2 - 2*b.se.m2, nrow=J, ncol=1)
upper <- as.matrix(b.hat.m2 + 2*b.se.m2, nrow=J, ncol=1)
b.hat <- as.matrix(b.hat.m2, nrow=J, ncol=1)
region.num.unorder <- as.matrix(regionvec, nrow=J, ncol=1)
region.nam.unorder <- as.matrix(uniq.name, nrow=J, ncol=1)

fig.mat <- as.data.frame(cbind(region.num.unorder, b.hat, lower, upper))
fig.mat <- as.data.frame(cbind(fig.mat, region.nam.unorder))

fig.mat <- fig.mat[order(fig.mat[,2]),]

region.num.order <- as.numeric(1:J)
region.num.order <- as.matrix(region.num.order, nrow=J, ncol=1)
fig.mat <- as.data.frame(cbind(fig.mat, region.num.order))


# Making the Figure
plot(fig.mat$region.num.order, fig.mat$V2, xlab= expression(paste("Region i, ranked by Estimated Treatment Effect ", beta[i])), ylab= expression(beta[i]), ylim=(c(-1.0,0.1)), xlim=(c(0,10)),xaxt='n')
#plot(fig.mat$region.num.order, fig.mat$V2, xlab= expression(paste("Region i, ranked by Estimated Treatment Effect ", beta[i])), ylab= expression(beta[i]))
segments(fig.mat$region.num.order, fig.mat$V3, fig.mat$region.num.order, fig.mat$V4, lwd=.5, col="gray10")
abline(h=0)
text(fig.mat$region.num.order, fig.mat$V2, fig.mat$region.nam.unorder, cex=1.3, pos=3, srt=45, offset = 0.48,  col="black") 
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_KStan_Drafts/regionspeccoeff_wctrl_logit.eps", horizontal=FALSE)
#	this is the one used in the MS



# Replicating the MLM with the ICC outcome measure

m2region_icc <- glmer(y_icc ~ x + (1+x|region) + age_under50 + male + postsec_educ + employed + income_av, family = binomial)

m2region_icc
ranef(m2region_icc)
se.ranef(m2region_icc)

model2_icc<-mtable(m2region_icc)
toLatex(model2_icc)

### Figure XX: Region specific treatment effects

b.hat.m2 <- coef(m2region_icc)$region[,2]
b.se.m2 <- se.coef(m2region_icc)$region[,2]
J <- 9
regionvec <- as.numeric(1:J)

# Make vector of country names
region.name <- as.vector(data$I1)
uniq.name <- sort(unique(region.name))

lower <- as.matrix(b.hat.m2 - 2*b.se.m2, nrow=J, ncol=1)
upper <- as.matrix(b.hat.m2 + 2*b.se.m2, nrow=J, ncol=1)
b.hat <- as.matrix(b.hat.m2, nrow=J, ncol=1)
region.num.unorder <- as.matrix(regionvec, nrow=J, ncol=1)
region.nam.unorder <- as.matrix(uniq.name, nrow=J, ncol=1)

fig.mat <- as.data.frame(cbind(region.num.unorder, b.hat, lower, upper))
fig.mat <- as.data.frame(cbind(fig.mat, region.nam.unorder))

fig.mat <- fig.mat[order(fig.mat[,2]),]

region.num.order <- as.numeric(1:J)
region.num.order <- as.matrix(region.num.order, nrow=J, ncol=1)
fig.mat <- as.data.frame(cbind(fig.mat, region.num.order))


# Making the Figure
plot(fig.mat$region.num.order, fig.mat$V2, xlab= expression(paste("Region i, ranked by Estimated Treatment Effect ", beta[i])), ylab= expression(beta[i]), ylim=(c(-0.6,0.5)), xlim=(c(0,10)),xaxt='n')
#plot(fig.mat$region.num.order, fig.mat$V2, xlab= expression(paste("Region i, ranked by Estimated Treatment Effect ", beta[i])), ylab= expression(beta[i]))
segments(fig.mat$region.num.order, fig.mat$V3, fig.mat$region.num.order, fig.mat$V4, lwd=.5, col="gray10")
abline(h=0)
text(fig.mat$region.num.order, fig.mat$V2, fig.mat$region.nam.unorder, cex=0.7, pos=3, offset = 0.48,  col="red") 
dev.print(postscript, "/Users/robertchaudoin/Dropbox/ICC Kyrgyzstan/CC_KStan_Drafts/regionspeccoeff_wctrl_logit_icc.eps", horizontal=FALSE)

